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Abstract. We study the transport of ultracold impurity atoms immersed in a Bose- 
Einstcin condensate (BEC) and trapped in a tight optical lattice. Within the strong- 
coupling regime, we derive an extended Hubbard model describing the dynamics of the 
impurities in terms of polarons, i.e. impurities dressed by a coherent state of Bogoliubov 
phonons. Using a generalized master equation based on this microscopic model we show 
that inelastic and dissipative phonon scattering results in (i) a crossover from coherent 
to incoherent transport of impurities with increasing BEC temperature and (ii) the 
emergence of a net atomic current across a tilted optical lattice. The dependence of 
the atomic current on the lattice tilt changes from ohmic conductance to negative 
differential conductance within an experimentally accessible parameter regime. This 
transition is accurately described by an Esaki-Tsu-type relation with the effective 
relaxation time of the impurities as a temperature-dependent parameter. 
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1. Introduction 

The study of impurities immersed in liquid helium in the 1960s opened a new chapter 
in the understanding of the structure and dynamics of a Bose- condensed fluid [1]. 
Two of the most prominent examples include dilute solutions of 3 He in superfluid 
4 He [1, 2] and the measurement of ionic mobilities in superfluid 4 He [1, 3]. More 
recently, the experimental realization of impurities in a Bose-Einstein condensate 
(BEC) [4,5] and the possibility to produce quantum degenerate atomic mixtures [6-11] 
have generated renewed interest in the physics of impurities. Of particular importance, 
collisionally induced transport of impurities in an ultra-cold bosonic bath has been 
recently observed [12] in an experimental setup of a similar type as to the one considered 
in this paper. 

In the context of liquid helium and Bose-Einstein condensates, a plethora of 
theoretical results have been obtained. Notably, the effective interaction between 
impurities [13-16] and the effective mass of impurities [17-21], both of which are induced 
by the condensate background, have been studied in detail. Moreover, the problem of 
self-trapping of static impurities has been addressed [22-24] in the framework of Gross- 
Pitaevskii (GP) theory [25,26], and quantized excitations around the ground state of 
the self-trapped impurity and the distorted BEC have been investigated in [27]. 



Figure 1. Impurities described by the localized Wannier functions Xj( r ) are trapped 
in the optical lattice potential Vt(r) and immersed in a nearly uniform BEC. A tilt 
imposed on the periodic potential leads to the formation of a Wannier-Stark ladder 
(dotted lines) with levels separated by the Bloch frequency wg. Inelastic scattering 
of Bogoliubov phonons results in a net atomic current across the tilted lattice. The 
current depends on the effective hopping J, the level separation lub and the effective 
relaxation time of the impurities. 



In the present paper we study the dynamics of impurities immersed in a nearly 
uniform BEC, which in addition are trapped in an optical lattice potential [28,29]. This 
setup has been considered in [30] by the present authors, however, we here motivate 
and develop the small polaron formalism in more detail and extend the results to the 
case of a tilted optical lattice potential, as illustrated in Fig. [TJ In contrast to similar 
previously-studied systems, the presence of an optical lattice allows us to completely 
control the kinetic energy of the impurities. Crucially, for impurities in the lowest Bloch 
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band of the optical lattice the kinetic energy is considerably reduced by increasing the 
depth of the lattice potential [28,29]. As a consequence, it is possible to access the 
so-called strong-coupling regime [31], where the interaction energy due to the coupling 
between the impurity and the BEC is much larger than the kinetic energy of the impurity. 
This regime is particularly interesting since it involves incoherent and dissipative multi- 
phonon processes which have a profound effect on the transport properties of the 
impurities. 

As a starting point we investigate the problem of static impurities based on 
the GP approximation and by means of a Bogoliubov description of the quantized 
excitations [26,32]. The main part of the paper is based on small polaron theory [33-35], 
where the polaron is composed of an impurity dressed by a coherent state of Bogoliubov 
phonons. This formalism accounts in a natural way for the effects of the BEC 
background on the impurities and allows us to describe the BEC deformation around 
the impurity in terms of displacement operators. As a principal result, we obtain 
an extended Hubbard model [36, 37] for the impurities, which includes the hopping 
of polarons and the effective impurity-impurity interaction mediated by the BEC. 
Moreover, the microscopic model is particularly suitable for the investigation of the 
transport of impurities from first principles within the framework of a generalized master 
equation (GME) approach [38-41]. 

The transport properties of the impurities are considered in two distinct setups. 
We first show that in a non-tilted optical lattice the Bogoliubov phonons induce a 
crossover from coherent to incoherent hopping as the BEC temperature increases [30]. 
In addition, we extend the results in [30] to the case of a tilted optical lattice (see Fig. [Q. 
We demonstrate that the inelastic phonon scattering responsible for the incoherent 
hopping of the impurities also provides the necessary relaxation process required for the 
emergence of a net atomic current across the lattice [42]. In particular, we find that 
the dependence of the current on the lattice tilt changes from ohmic conductance to 
negative differential conductance (NDC) for sufficiently low BEC temperatures. So far, 
to our knowledge, NDC has only been observed in a non-degenerate mixture of ultra-cold 
40 K and 87 Rb atoms [12] and in semiconductor superlattices [43,44], where the voltage- 
current dependence is known to obey the Esaki-Tsu relation [45]. By exploiting the 
analogy with a solid state system we show that the transition from ohmic conductance 
to NDC is accurately described by a similar relation, which allows us to determine the 
effective relaxation time of the impurities as a function of the BEC temperature. 

The implementation of our setup relies on recent experimental progress in the 
production of quantum degenerate atomic mixtures [6-9] and their confinement in 
optical lattice potentials [10, 11]. As proposed in [46], the impurities in the optical 
lattice can be cooled to extremely low temperatures by exploiting the collisional 
interaction with the surrounding BEC. The tilt imposed on the lattice potential may be 
achieved either by accelerating the lattice [47] or superimposing an additional harmonic 
potential [12,48]. 
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An essential requirement for our model is that neither interactions with impurities 
nor the trapping potential confining the impurities impairs the ability of the BEC 
to sustain phonon-like excitations. The first condition limits the number of impurity 
atoms [10, 11,49], and thus we assume that the filling factor of the lattice is much lower 
than one, whereas the second requirement can be met by using a species-specific optical 
lattice potential [50]. Unlike in the case of self-trapped impurities [22-24], we assume 
that the one-particle states of the impurities are not modified by the BEC, which can 



be achieved by sufficiently tight impurity trapping as shown in Appendix C 

The paper is organized as follows. In Section 2 we present our model. In Section 3, 
we motivate the small polaron formalism by considering the problem of static impurities 
based on the GP equation and within the framework of Bogoliubov theory. In Section 4, 
the full dynamics of the impurities is included. We develop the small polaron formalism 
in detail and derive an effective Hamiltonian for hopping impurities. In Section 5, 
the GME is used to investigate the transport properties of the impurity atoms. First, 
we discuss the crossover from coherent to incoherent hopping and subsequently study 
the dependence of a net atomic current across a tilted optical lattice on the system 
parameters. We discuss possible extensions of our model and conclude in Section 6. 
Throughout the paper we consider impurity hopping only along one single direction for 
notational convenience. However, the generalization to hopping along more than one 
direction is straightforward. 

2. Model 

The Hamiltonian of the system is composed of three parts H = Hb + H x + Hj, where 
Hb is the Hamiltonian of the Bose-gas, H x governs the dynamics of the impurity atoms 
and Hi describes the interactions between the impurities and the condensate atoms. 

The Bose-gas is assumed to be at a temperature T well below the critical 
temperature T c , so that most of the atoms are in the Bose-condensed state. The 
boson-boson interaction is represented by a pseudo-potential g5(r — r'), where the 
coupling constant g > depends on the s-wave scattering length of the condensate 
atoms. Accordingly, the grand canonical Hamiltonian Hb of the Bose-gas reads 

Hb — J dr^t(r) [h - fi b ] $(r) + || dr^(r)^( r )^(r)^(r) , (1) 

where ^ is the chemical potential and 

H = -^-V 2 + V ext (r) (2) 

is the single-particle Hamiltonian of the non-interacting gas, with m;, the mass of a boson 
and V r cxt (r) is a shallow external trapping potential. The boson field operators $(r) and 
ip(r) create and annihilate an atom at position r, respectively, and satisfy the bosonic 
commutation relations [V>(r), V^(r')] = 5(r — r') and [^(r), -0(r')] = [^(r),^^')] = 0. 
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To be specific we consider bosonic impurity atoms with mass m a loaded into an 
optical lattice, which are described by the impurity field operator x(r). The optical 
lattice potential for the impurities is of the form Vl(t) = V x sm 2 (kx) + V y sin 2 (ky) + 
V z sin 2 (kz), with the wave- vector k = 2tt/X and A the wavelength of the laser beams. 
The depth of the lattice V% (£ = x, y, z) is determined by the intensity of the 
corresponding laser beam and conveniently measured in units of the recoil energy 
E R = h 2 k 2 /2m a . Since we are interested in the strong-coupling regime we assume that 
the depth of the lattice exceeds several recoil energies. As a consequence, the impurity 
dynamics is accurately described in the tight-binding approximation by the Bose- 
Hubbard model [28,29], where the mode-functions of the impurities Xj{ r ) are Wannier 
functions of the lowest Bloch band localized at site j, satisfying the normalization 
condition / dr|xj(r)| 2 = 1. Accordingly, the impurity field operator is expanded as 
X.( r ) = J2j Xj{ r )®j an d the Bose-Hubbard Hamiltonian is given by 

H x = + \ U Y^ %(% ~ 1 )+l J 'aJ2 % + huJ B^Zj flj , (3) 

3 3 3 

where \i a is the energy offset, U is the on-site interaction strength, J is the hopping 
matrix element between adjacent sites and denotes the sum over nearest 

neighbours. The operators a] (a,j) create (annihilate) an impurity at lattice site j and 
satisfy the bosonic commutation relations [cij,a]] = Sij and [a*, a,] = = 0, and 

hj = alaj is the number operator. The last term in Eq. ([3]) was added to allow for 
a tilted lattice with the energy levels of the Wannier states Xj{ r ) separated by Bloch 
frequency ub, as shown in Fig. (pQ). 

For the lattice depths considered in this paper it is possible to expand the potential 
wells about their minima and approximate the Wannier functions by harmonic oscillator 
ground states. The oscillation frequencies of the wells are u>£ = 2V EnVi/h and 
the spread of the mode- functions Xj{ r ) is given by the harmonic oscillator length 
at = \Juj m a uj£ . Explicitly, the Gaussian mode-function in one dimension reads 

X 3 Ax) = (7n7 2 )- 1/4 exp [-{x - x J ) 2 /{2a 2 )] , (4) 

where Xj = aj with a = A/2 being the lattice spacing. We note that the mode- functions 
Xj,<r( r ) are highly localized on a scale short compared to the lattice spacing a since 
aij a = (Vi/ ' Er)~ 1 ^ /tt 1 for sufficiently deep lattices. Consistency of the tight-binging 
model requires that excitations into the first excited Bloch band due to Landau-Zener 
tunneling or on-site interactions are negligible, i.e. ojb <C uj£ and 1/2 Urijirij — 1) <C fru>e, 
with rij the occupation numbers. These inequalities are readily satisfied in practice. 

The density-density interaction between the impurities and the Bose-gas atoms is 
again described by a pseudo-potential k8(t — r'), with k the coupling constant, and the 
interaction Hamiltonian is of the form 

Hi = kJ drxt( r )x( r )^t( r )^( r ). (5) 
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We note that the effect of Hb and Hi dominate the dynamics of the system in the strong- 
coupling regime and hence the kinetic part in H x can be treated as a perturbation. 



3. Static impurities 

In a first step we investigate the interaction of static impurities with the BEC, i.e. we 
neglect the hopping term in H x . We first derive a mean-field Hamiltonian based on the 
GP approximation and subsequently quantize the small-amplitude oscillations about 
the GP ground state using standard Bogoliubov theory [26,32]. 

In the GP approximation, the bosonic field operator ip(r) is replaced by the classical 
order parameter of the condensate ip(r). Accordingly, the impurities are described by 
the impurity density p x (r) = (j^ (r)x(r)) = J2j n j\Xj( r )\ 2 , where the average is taken 
over a fixed product state |T) representing a set of static impurities, i.e. |T) = aj|0), 
with {j} a set of occupied sites and |0) the impurity vacuum. We linearize the equation 
for ip(r) by considering small deviations i?(r) = ip(r) — ipo(r) from the ground state of 
the condensate ipo(r) in absence of impurities. This approach is valid provided that 
|i?(r)|/i/)o( r ) ^ 1 and corresponds to an expansion of Hb + Hi to second order in k. 
Explicitly, by replacing ip(r) with ^>o( r )+$( r ) i n Hb + Hi we obtain the GP Hamiltonian 



Hop 



H^ + H# + Hun with 



drk*(r)^(r) 



MMr)\ 2 + |l^o(r)| 4 \ + k I drp x (r)|^o(r)| 2 



(6a) 



H*= Jdr {tf*(r) [h - p b + 2^o(r)| 2 ] #(')} 

+ 9 - [ dr{r(r)[^ (r)] 2 r(r) + ^(r)[^(r)]^(r)} , 



(66) 



F lin = K J dr p x (r) [Vo(r)r (r) + VoW^W] . (6c) 

In absence of impurities, i.e. for k = and i9(r) identically zero, the stationary condition 
5H GP / 5iPq{y) = implies that ^o( r ) satisfies the time-independent GP equation 

^o + ^|^o(r)| 2 ]^o(r)=^ (r). (7) 

The deformed ground state of Hqp due to the presence of impurities is found by imposing 
the stationary condition SHap/S"d*(r) = or equivalently 

'H -p b + 2^ (r)| 2 ] tf(r) + gM*)?#* W + *P x (r)Mr) = , (8) 

which is the Bogoliubov-de Gennes equation [26,32] for the zero-frequency mode. 

For the case of a homogeneous BEC, the effect of the impurities on the ground 
state of the condensate can be expressed in terms of Green's functions [22,23,27]. If 
we assume for simplicity that both ^o( r ) and $(r) are real then Eq. (jSj) reduces to the 
modified Helmholtz equation 

21 



0(r) 



2k 



DPx( r ) 



(9) 
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Figure 2. The order parameter ip(r) of a two-dimensional homogeneous condensate 
deformed by two highly localized impurities with n > according to Eq. (fT0|) . The 



deformation of the BEC results in an effective interaction potential Vij between the 
impurities at lattice sites i and j. The range of the potential is characterized by the 
healing length £, which is comparable to the lattice spacing a for realistic experimental 
parameters. 



where £ = h/ ^Jm^gn^ is the healing length, D is the number of spatial dimensions, and 
n o = iV'ol 2 is the condensate density. The solution of Eq. (EJ) is found to be [51] 



0(r) 



K 



where the Green's functions are 

1 



J dr'3(r-r> x (r' 



£m(r) 
^2D(r) 
^ 3D (r) 



2 exp(-2|r|/£) , 

^K (2|r|/0, 
1 exp(-2|r|/0 



(10) 



'11a) 



[lib) 



[11c) 



27T \v\/i 

with K Q (x) the modified Bessel function of the second kind. Thus, independent of the 
dimension D of the system, the BEC deformation induced by the impurities falls off 
exponentially on a length scale set by £, as illustrated in Fig. (T5]). 

It follows from Eq. ([TO]) that for occupation numbers rij ~ 1 the condition 
|$(r)|/?/> (r) <c 1 is equivalent to 

D 

( .= — |-| <1, (12) 




with d 



71, 



-l/D 



the average separation of the condensate atoms. The generalization of 



this condition to the non-homogeneous case is a(r) = (|ft|/<?) [c?(r)/£(r)] D <^ 1 under 
the assumption that the BEC is nearly uniform. We note that a(r) oc no(r)^ D_2 ^ 2 , 
and hence a(r) diverges in ID as rio( r ) ~~ * 0, e.g. near the boundary of the condensate. 
However, this is consistent with the GP approach, which is no longer applicable in the 
dilute limit f (r)/d(r) < 1 of a ID Bose-gas [26]. 
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The change in energy of the BEC due to the impurities is found by inserting the 
formal solution for $(r) in Eq. (1101) into H G p- Provided that $(r) satisfies Eq. (jSJ) the 
identity H\ in + 2H$ = 0, or equivalently H GP = H^ + ^H\ in , holds. Using the latter 
expression for Hqp we find 

H G p = J2 - e i) n i - £ f i n i ( n i ~ X ) ~ 2 ^ Vi,i ^ ' ^ 

where £ = /m is the first order contribution to the mean-field shift. In Eq. ([TBI 
we have neglected constant terms which do not depend on the impurity configuration, 
i.e. contributions containing tpo(r) only. The off-site interaction potential between the 
impurities is given by 

Vy = ^J drdr'| Xi (r)| 2 ^(r-r')|x,(r')| 2 , (14) 

and £j = \Vjj is the potential energy of an impurity, both resulting from the deformation 
of the condensate, as illustrated in Fig. (J2J). 

In order to quantize the GP solution we consider small excitations £(r) of the 
system, which obey bosonic commutation relations, around the static GP ground state 
ipo(r) +i?(r) of the condensate [26,32]. This corresponds to an expansion of the bosonic 
field operator as ip(r) = ^o( r )+^( r )+C( r )- The Hamiltonian governing the evolution 
of C(r) can be obtained by substituting #(r) + £(r) for i?(r) in Hgp- By collecting the 
terms containing £(r) and C^(r) we find 

H c = Jdr {C f (r) [Ho - & + 2<#o(r)| 2 ] f (r)} 

+ f /dr {ft(p)[^(r)] a Ct(r)+ C(r)[V>o*(r)] 2 C(r)} , (15) 

where the linear terms in £(r) and C^( r ) vanish identically since $(r) satisfies Eq. (jHJ). 
The Hamiltonian in Eq. (fl~5l) can be diagonalized by the standard Bogoliubov 
transformation [26,32] 

C(r) = £[Mr)/i + <(r)/3t] . (16) 

V 

Here, u v (r) and v*(r) are complex functions and the operators pi (/3„) create (annihilate) 
a Bogoliubov quasi-particle, with quantum numbers v, and satisfy bosonic commutation 
relations. The Bogoliubov transformation in Eq. (1161) reduces the Hamiltonian H^ to a 
collection of noninteracting quasi-particles provided that the coefficients u v (r) and v v (r) 
obey the Bogoliubov-de Gennes equations [26,32] 



H - p, b + 2g\ijj (r)\ 2 u„(r) + g[%{r)fv v (r) = hu u u u 



H Q - fi b + 2g\i[> (r)\ v u (r) + g[ipo(r)) u v (r) = -Tiuj u v 



(17) 



with the spectrum tku v . Under the condition that Eqs. ffTTj) are satisfied one finds that 
H^ takes the form H^ = yZ u Ku u /3lj3 u , where constant terms depending on ipo(r) and 
v v {t) only were neglected [j. 

X The details of establishing the above form of are contained in Ref. [32] . 
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Consequently, the total Hamiltonian H stat of the BEC and the static impurities is 
composed of three parts 

-Hstat = H{ + H GP + (H x ) , (18) 

where Hq governs the dynamics of the Bogoliubov quasi-particles and Hop is the GP 
ground state energy, which for the case of a homogeneous BEC takes the simple form 
in Eq. 0131) . The third term (H x ) represents the average value of H x with respect to 
the fixed product state |T) introduced earlier with rij impurities at each site j, giving 
explicitly 

( H x) = \ u Y. n i ( n j - l )+^aY. n i + huJ B 3 %' • (19) 

3 3 3 

The ground state of the system corresponds to the Bogoliubov vacuum defined by 
/3^|vac) = 0. 



4. Hopping impurities in the polaron picture 



Given the Hamiltonian H sta x, derived in the previous section it is straightforward to 
determine the ground state energy for a set of static impurities. However, an analysis 
of the full dynamics of the impurities requires an alternative description in terms of 
polarons, i.e. impurities dressed by a coherent state of Bogoliubov quasi-particles. In this 
picture small polaron theory allows us, for example, to calculate the effective hopping 
matrix element for the impurities, which takes the BEC background into account. 

Our approach is based on the observation that the Hamiltonian of the quantized 
excitations and consequently the Bogoliubov-de Gennes equations (fT7j) are independent 
of p x (r). In other words, the effect of the impurities is only to shift the equilibrium 
position of the modes j3 v without changing the spectrum hui u . Therefore it is possible to 
first expand the bosonic field operator as $(r) = ^o( r ) + ^( r )> with $(r) = $(r) + C( r )> 
subsequently express i?(r) in terms of Bogoliubov modes about the state ipo( T ) an d 
finally shift their equilibrium positions in order to (approximately) minimize the total 
energy of the system. 

Similarly to the case of static impurities, we first replace ?/>(r) with ^>o( r ) + "^( r ) m 
Hp + Hi to obtain the Hamiltonian Hqp = H^ + H$ + Hn n with 

"dr^WxWI^oWl 2 , (20a) 



H, 



li'O 



dr^oW^o(r) -fi b \M*)\' 



fl^oWl 4 



+ re 



%= f drp( r ) [H -fi b + 2g\Mr)\ 2 ]Hr)} 



+ ||drp(r)[^(r)] a ^(r)+^(r)[V€(r)]^(r)} , 



Fiin = k / drx f (r)x(r) Uo(r)i? f (r) +^o( r )^( r 



(206) 
(20 c) 
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We note that now H G p contains the density operator xK r )x( r ) instead of p x (r) in order 
to take into account the full dynamics of the impurities. The expansion of $(r) in terms 
of Bogoliubov modes reads 

$(r) =Y[u v (r)K+v:(r)bl\ , (21) 

V 

where the spectrum hu u and coefficients u v {r) and v u (r) are determined by Eqs. ffT7|) . 
The bosonic operators 6J, (b u ) create (annihilate) a Bogoliubov excitation around the 
ground state of the condensate ipo(r) in absence of impurities, and thus, importantly, 
do not annihilate the vacuum |vac) defined in the previous section, i.e. 6j,|vac) 7^ 0. By 
substituting the expansion in Eq. (|2"T|) for $(r) in the Hamiltonian Hgp and using the 
identity xH r )x( r ) = J2i,j Xi ( r )Xj( r ) a I a j we find that up to constant terms T 



<l'o 



E4 



'■J 



H# = J2huj u blK, 

V 

H hn = E hu u [M i;ji X + M* jjU bl 



with the matrix elements 

= K 

K 



£■ ■ 



(22 a) 
{22b) 
(22 c) 

(23) 

where we assumed for simplicity that V ; o( r ) is rea l- The non-local couplings M^ v 
and Si j with i ^ j resulting from the off-diagonal elements in x ( r )x( r ) are highly 
suppressed because the product of two mode-functions x«( r ) an d Xj( r ) with i ^ j is 
exponentially small. As a consequence, the evolution of the BEC and the impurities, 
including the full impurity Hamiltonian H x , is accurately described by the Hubbard- 
Holstein Hamiltonian [33-35] 



drn (r)x*( r )Xj( r ) , 
/ dr^ (r)[^(r) + ^(r)]x*( r )Xj( r ) , 



3>v 



M, 



where we discarded the non-local couplings and introduced Bj = £jj and Mj tV 

Since we intend to treat the dynamics of the impurities as a perturbation we shift 
the operators b v and b\ in such a way that H ho \ in Eq. (|24|) is diagonal, i.e. the total 
energy of the system is exactly minimized, for the case J = 0. This is achieved by the 
unitary Lang-Firsov transformation [34, 35] 



U = exp 



E ( M lM - M h X 



n : 



3V 



(25) 



from which it follows that by applying the Baker- Campbell-Hausdorff formula [52] 



Ubllfl = bt - M, hV h h UhjW 



n. 



and Uatfr = aW, with 



exp 



E (M*Jt - M h X 



(26) 



As for , the details of establishing the above form of H$ are contained in Ref. [32] . 
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The operator Jtj is the sought after displacement operator that creates a coherent state 
of Bogoliubov quasi-particles, i.e. a condensate deformation around the impurity. As a 
result, the transformed Hamiltonian Hlf = UH^iU^ is given by 

h lf = - Jj2(XAy(Xj°<j) + ^ E % " X ) 

+ E P-i % + huJ B E i ™i _ g E + E huJ M^ > ( 27 ) 

with the effective energy offset /2j = /i a + £j- — Ej, the effective on-site interaction 
strength Uj = U — 2Ej, the interaction potential V itj = Y,u^v {M iiV M* y + M* u Mj jV J, 
and Ej = \Vjj the so-called polaronic level shift [34,35]. 

In the case of static impurities and, more generally, in the limit ( = J/£j < 1, 
the polarons created by olxj are the appropriate quasi-particles. Consequently, the 
Hamiltonian in Eq. ([27)1 describes the dynamics of hopping polarons according to an 
extended Hubbard model [36,37] with an non-retarded interaction potential Vij. The 
contributions to Hlf are qualitatively the same as for static impurities, except for the 
additional hopping term. In the thermodynamic limit, the interaction potential Vij and 
the polaronic level shift Ej are identical, respectively, to their GP counterparts Vy and 
Ej. The reason for the simple connection between GP theory and small polaron results 
is that both involve a linearization of the equations describing the condensate. It should 
be noted that the interaction potential Vy and the polaronic level shift Ej can also 
be obtained exactly from H hol by applying standard Rayleigh-Schrodinger perturbation 
theory up to second order in k since all higher order terms vanish. However, the merit 
of using the Lang-Firsov transformation lies in the fact that it yields a true many-body 
description of the state of the system, which would require a summation of perturbation 
terms to all orders in k [53]. 

We gain qualitative and quantitative insight into the dependence of the quantities 
in Eq. fT27j) on the system parameters by considering the specific case of a homogeneous 
BEC with the total Hamiltonian 

Hlf = - J (*A) + luJ2hj(hj - 1) 

+ pL E flj + huJ B E j ''TLj - ~ E Vi,j hiKj + E ^q&q&q , ( 28 ) 

with the energy offset = + kti — E p , the on-site interaction strength U = U — 2E P , 
the polaronic level shift Ej = E p and the phonon momentum q. For a homogeneous BEC 
the Bogoliubov coefficients are of the form u q (r) = w q exp(iq-r) and f q (r) = w q exp(iq-r) 
with [26,32] 

Here, e q = (hq) 2 /2mb is the free particle energy, hujq = we q (e q + 2(?no) is the 
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Figure 3. (a) The function — G^jj(r, a) plotted versus r for er/£ = 0.1 (solid line), 
cr/£ = 0.15 (dashed line), a/£ = 0.2 (dotted line). The potential V(r) oc G 3D (r,a) 
falls off on a scale set by £ and has an increasing depth V(0) with decreasing a. (b) 
The function — G(0,a) plotted versus a for 3D (solid line), 2D (dashed line) and ID 
(dotted line). The depth of the potential V(0) oc G(0, a) and £ p — \V{0) depend 
strongly on a in 3D and only weakly in 2D and ID. The relevant range of a can be 
estimated for deep lattices by er/£ ~ (a/ '£) /[n(Ve/ 'Er) 1 ^]. 



Bogoliubov dispersion relation and VL is the quantization volume. The corresponding 
matrix elements are found to be 



Mj < q= V(^ /j(q) ' (30) 

where 

/i(q) = -j=J dr|xj(r)| 2 exp(iq- r) . (31) 
For the Gaussian mode-function Xj,<r( r ) m Eq. (HI) the factor /j(q) takes the simple form 

/i(q) = ~m II ex P (-9j°i/ 4 ) ex P (iq ■ t j) ■ ( 32 ) 



n 

The q-dependence of the impurity-phonon coupling is Mj q oc fj (q) / y^qj in the long 
wavelength limit |q| <C l/£, which corresponds to a coupling to acoustic phonons 
in a solid state system [35], whereas in the free particle regime |q| 3> l/£ we have 

The potential Vij and the polaronic level shift E p reduce to a sum over all 
momenta q and can be evaluated in the thermodynamic limit J2q ~^ (2vr) _j:, /dq 
where V!y — > V{j. In particular, for the Gaussian mode-functions Xj,a( r ) one finds 

Hi^^J^-r^), (33) 

where the functions G(r, a) are defined in Appendix A| and plotted in Fig. [3^ for a 
three-dimensional system. It follows from the definition of G(r, 0) that Q(r) = G(r, 0), 
and hence the shape of the potential is determined by the Green's function Q(r) in 
the limit u € (. The cr-dependence of G(0,a), which determines the depth of the 
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interaction potential Vij and the level shift £ p = ^Vj t j, is shown in Fig. [3b. The range 
of the potential, characterized by the healing length £, is comparable to the lattice 
spacing a for realistic experimental parameters, and hence the off-site terms Vjj+i are 
non- negligible. In particular, as shown in [30,54], the off-site interactions can lead to 
the aggregation of polarons on adjacent lattice sites into stable clusters, which are not 
prone to loss from three-body inelastic collisions. 



5. Transport 

The coupling to the Bogoliubov phonons via the operators Xj and Xj changes the 
transport properties of the impurities notably. Since we assume the filling factor of 
the lattice to be much lower than one we can investigate the transport properties by 
considering a single polaron with the Hamiltonian ho + hi given by 

h = hu B 3 + huJ AK ' ( 34 ) 
i q 

h I = -J^2(X l a i )\X J a J ). (35) 

(ij) 

To start we investigate the crossover from coherent to diffusive hopping [33, 35] in a 
non-tilted lattice {oj b = 0) and then extend the result to a tilted lattice (u B ^ 0) to 
demonstrate the emergence of a net atomic current across the lattice [42] due to energy 
dissipation into the BEC. 



5.1. Coherent versus incoherent transport 

We first consider coherent hopping of polarons at small BEC temperatures k B T -C E p , 
where incoherent phonon scattering is highly suppressed. In the strong-coupling regime 
( = J/Ej -C 1, and hence the hopping term in Eq. (155]) can be treated as a perturbation. 
The degeneracy of the Wannier states requires a change into the Bloch basis 

l k > = 7^£ ex P( ik - r i)«]l >> (36) 
ViV j 

with k the quasi-momentum and N the number of lattice sites. Applying standard 
perturbation theory in the Bloch basis and describing the state of the system by 
|k, {iV q }), where {iV q } is the phonon configuration with phonon occupation numbers 
iV q , we find the polaron energy up to first order in £ 

E{k) = jl a - 2Jcos(k-a) , (37) 

where we defined the effective hopping J = JX) q (-^ r q|AjA > :;+ i|A^ q ) and a is the position 
vector connecting two nearest neighbor sites. In particular, for the case of a thermal 
phonon distribution with occupation numbers iV q = (e h ^il k B T _ x) _1 the effective 
hopping is [33-35] 

J = J exp { - J2 l M o, q | 2 [1 - cos(q-a)] (2iV q + 1) } . (38) 
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Thus, the hopping bandwidth of the polaron band is highly suppressed with increasing 
coupling k and temperature T. 

At high temperatures E p <C k B T <C ksT c inelastic scattering, in which phonons are 
emitted and absorbed, becomes dominant, and thus the transport of impurities through 
the lattice changes from being purely coherent to incoherent. While matrix elements 
({Nq}\XjXj + x\{Nqy) involving two different phonon configurations {iV q } and {Nq}' 
vanish at zero temperature they can take non-zero values for T > 0. The condition 
for energy and momentum conservation during a hopping event implies that incoherent 
hopping is dominated by a three-phonon process that involves phonons with a linear 



dispersion Tluj^ = hc\q\, where c ~ ygno/rrib is the speed of sound. This process is 
reminiscent of the well-known Beliaev decay of phonons [26] . 

We investigate the incoherent transport properties by using a generalized master 
equation (GME) [38] for the site occupation probabilities Pj(t) of the impurity. This 
formalism has been applied to the transfer of excitons in the presence of electron-phonon 
coupling [39,40] and is based on the Nakajima-Zwanzig projection method [41]. The 
generalized master equation is of the form 

= [ ds HWiAs)lPi(t -s)- m - a)] , (39) 

where the effect of the condensate is encoded in the memory functions Wij(s). As shown 



in 



Appendix B the memory function to second order in ( is given by 

2 



Wij(s) = 28 j>i±1 - Re 



exp{-2]T|Mo, q | 2 [l-cos(q-a)] 



q^O 



x [(JV q + 1)(1 - e^ s ) + iV q (l - e^)]} exp(±k^s) 



(40) 



In the case ujb = 0, the nontrivial part of Wi t j(s) takes the values 2(J/h) 2 at s = and 
2(J/h) 2 in the limit s — > oo due to the cancellation of highly oscillating terms, as shown 
in Fig. UK 

In the regime k%T <ti E p , the effective hopping J is comparable to J and the 
memory function Wij is well approximated by 2(J/h) 2 Q(s), with O(s) the Heaviside step 
function, and thus describes purely coherent hopping. In the regime E p -C k B T /cbT c , 
processes involving thermal phonons become dominant and coherent hopping is highly 
suppressed, i.e. J <^ J. In this case the memory function W it j drops off sufficiently fast 
for the Markov approximation to be valid, as illustrated in Fig. S^. More precisely, one 
can replace Pj(t — s) by Pj(t) in Eq. fl3"9l and after intergration over s the GME reduces 
to the standard Pauli master equation 
dPi(t) 



dt j 



X>«ft(*)-^(*)]> (41) 



where the hopping rates Wi t j are given by 

/ dsiW^s) - limWiM ■ (42) 
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Figure 4. Coherent and diffusive hopping in a one-dimensional system: (a) The 
memory function h Wij(t)/(2J 2 ) plotted versus time for fc^T = (dotted line), 
ksT = 5E p (dashed line) and k B T = lbE p (solid line). The memory function drops 
off rapidly for k B T 3> E p , indicating the dominance of incoherent hopping, (b) - 
(d) The evolution of the occupation probabilities Pj (t) of a single impurity initially 
localized at site j = according to the GME with the memory functions in (a), (b) For 
small temperatures k B T <C E p the hopping is coherent with two wave-packets moving 
away from j = 0. (d) For high temperatures k B T 3> E p the inelastic scattering 
of phonons results in a diffusive motion of the impurity, where the probability Pj(t) 
remains peaked at j = 0. The lattice with spacing a = 395nm and J = 2.45 x 10~ 2 Er 
contains a single 41 K atom. The BEC with d = 200nm and £ = 652nm consists of 
87 Rb atoms; n/g = 2.58, E p /k B » lOnK and J/% w 1.2kHz. 



The Pauli master equation describes purely incoherent hopping with a thermally 
activated hopping rate Wij [33,35]. 

The evolution of an initially localized impurity at different temperatures T for a 
one-dimensional 41 K- 87 Rb system [8] is shown in Figs. Hb -HH, which were obtained 
by numerically solving the GME with the memory function Wij(s) in Eq. (j40l) . It can 
be seen that for small temperatures ksT E p the hopping is coherent with two wave- 
packets moving away from the initial position of the impurity atom. In contrast, for high 
temperatures k^T ^> E p inelastic scattering of phonons results in a diffusive motion of 
the impurity, where the probability Pj(t) remains peaked at the initial position of the 
impurity. 




Figure 5. Analysis of the mean-square displacement l 2 (t) = J^i l 2 Pi(t)i obtained from 
the evolution of an initially localized impurity in a one-dimensional system for the time 
tcvo\ = Wti/J according to the GME. The mean-squared displacement was assumed 
to be of the form l 2 (t) — At a . (a) The exponent a versus temperature T for the full 
evolution time (solid line) and the period t evo i/2 to t evo \ (dashed line). The drop from 
a » 2 at zero temperature to a ~ 1 at high temperatures kgT 3> E p clearly indicates 
the crossover from coherent to diffusive transport, (b) The prefactor A in units of 
(J/h) a versus temperature T for the full evolution time (solid line) and the period 
tevoi/2 to i ovo i (dashed line). The system parameters are the same as for Fig. |4] 



This crossover from coherent to diffusive hopping can be quantitatively analyzed 
by considering the mean-squared displacement of the impurity, l 2 (t) = J2i l 2 Pi(f)i which 
we assume to be of the form l 2 {t) = At a . The exponent a takes the value a = 2 for 
a purely coherent process, whereas a = 1 for diffusive hopping. Figure ([SD shows a 
and A in units of (J/h) a as functions of the BEC temperature T, which were obtained 
from the evolution (according to the GME) of an impurity initially localized at j = 0. 
We see that the exponent a drops from a ~ 2 at zero temperature to a ~ 1 at high 
temperatures ksT 3> E p , thereby clearly indicating the crossover from coherent to 
diffusive transport. We note that the transition from coherent to diffusive hopping 
takes place in a temperature regime accessible to experimental study and therefore, 
importantly, may be observable. 

5.2. Atomic current across a tilted lattice 

The inelastic phonon scattering responsible for the incoherent hopping of the impurities 
also provides the necessary relaxation process required for the emergence of a net atomic 
current across a tilted optical lattice. This is in contrast to coherent Bloch oscillations, 
which occur in an optical lattice system in absence of incoherent relaxation effects 
or dephasing [47]. As pointed out in [42], the dependence of the atomic current on 
the lattice tilt Hub changes from ohmic conductance to NDC in agreement with the 
theoretical model for electron transport introduced by Esaki and Tsu [45]. 

To demonstrate the emergence of a net atomic current, and, in particular, to show 
that impurities exhibit NDC, we consider the evolution of a localized impurity atom 
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Figure 6. Current-voltage relation for a one-dimensional system: (a) The drift velocity 
Vd in units of vq = Ja/% as a function of the lattice tilt Hlob/J for temperatures 
k B T = (+), k B T = 5E p (o) and k B T = 15E P (x) according to the GME and the 
best fit of the Esaki-Tsu-type relation in Eq. (|43[) (dotted lines) . The dependence of the 
current on the tilt changes from ohmic conductance Vd ~ Tiujb to negative differential 
conductance Vd ~ 1/(Hu>b)- (b) The relaxation time r (o) in units of To = h/griQ 
and the prefactor 7 (+) yielding the best fit of the Esaki-Tsu-type relation to the 
numerical results. The relaxation time r of the impurities decreases with increasing 
BEC temperature T and the prefactor 7 varies only slightly since the exponential 
temperature dependence of the current is accounted for by vq. The system parameters 
are the same as for Fig. [4] 



in a one-dimensional system. With the impurity initially at site j = we determine 
its average position Xd = J2j aj Pj(td) after a fixed drift time td of the order of K/J. 
This allows us to determine the drift velocity Vd = Xdftd as a function of the lattice tilt 
Thub and the temperature T of the BEC. In analogy with a solid state system, the drift 
velocity Vd and the lattice tilt Tiuib correspond to the current and voltage, respectively 

Figure [6^ shows the voltage-current relation at different temperatures T, which was 
obtained by numerically solving the GME with the memory function in Eq. fHUl) . We 
see that for a small lattice tilt the system exhibits ohmic behavior Vd ~ fruiB, whereas 
for a large lattice tilt the current decreases with increasing voltage as Vd ~ l/(hu B ), 
i.e. the impurities feature NDC. 

Following [42] we describe the voltage-current relation for the impurities by an 
Esaki-Tsu-type relation 

^ = 2 ™TT& (43) 

where vq = Ja/% the characteristic drift velocity, r the effective relaxation time of 
the impurities and 7 a dimensionless prefactor. Fitting the Esaki-Tsu-type relation 
in Eq. (|43l) to the numerical data allows us to extract the parameters r and 7, both 
depending on the BEC temperature. As can be seen in Fig. [6b, the effective relaxation 
time t decreases with increasing BEC temperature T and is of the order of To = h/(gno). 
The significance of 1/r is that of an average collision rate between the impurities and 
Bogoliubov excitations, which would allow us, for example, to formulate the problem of 
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transport in terms of a classical Boltzmann equation for the distribution function of the 
impurities [35]. The prefactor 7 varies only slightly since the exponential dependence 
of the current on the temperature is accounted for by Vq. We note that independently 
of ujbt the maximum drift velocity is given by 7^0- 

The Esaki-Tsu-type relation in Eq. (1431) reflects the competition between coherent 
and incoherent, dissipative processes. In the collision dominated regime ujbt "C 1, 
inelastic scattering with phonons destroys Bloch oscillations, whereas in the collisionless 
regime ubt ^> 1 the evolution of the impurities is mainly coherent, i.e. Bloch oscillations 
of the impurities lead to a suppression of the net current. The crossover between the two 
regimes is most pronounced at zero temperature, where J is comparable to J. However, 
the change from ohmic to negative differential conductance is identifiable even at finite 
temperatures and thus should be observable in an experimental setup similar to the one 
used in [12]. 

6. Conclusion 

We have studied the transport of impurity atoms in the strong-coupling regime, where 
the interaction energy due to the coupling between the BEC and the impurities 
dominates their dynamics. Within this regime, we have formulated an extended 
Hubbard model describing the impurities in terms of polarons, i.e. impurities dressed by 
a coherent state of Bogoliubov phonons. The model accommodates hopping of polarons 
and the effective off-site impurity-impurity interaction mediated by the BEC. 

Based on the extended Hubbard model we have shown from first principles that 
inelastic phonon scattering results in a crossover from coherent to incoherent hopping 
and leads to the emergence of a net atomic current across a tilted optical lattice. In 
particular, we have found that the dependence of the current on the lattice tilt changes 
from ohmic conductance to negative differential conductance for sufficiently low BEC 
temperatures. Notably, this transition is accurately described by an Esaki-Tsu-type 
relation with the effective relaxation time of the impurities as a temperature-dependent 
parameter. 

Using the techniques introduced in this paper, qualitatively similar phenomena 
can also be shown to occur for fermionic impurities and, moreover, for impurities of 
different species [55]. For instance, in the case of two impurity species A and B with 
the couplings ka > and Kb < 0, respectively, the effective off-site impurity-impurity 
interaction is attractive for the same species, but repulsive for different species. In 
either case, observation of the phenomena reported in this paper lies within the reach of 
current experiments, which may give new insight into the interplay between coherent, 
incoherent and dissipative processes in many-body systems. 
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Appendix A. Definition of the functions G(r, a) 

The function G(r, a), which are a generalization of the Green's functions Q(r) for o > 0, 
are defined by 

2" 



Cm r, a) = - / dy exp 

7r J-oo y 2 + 4 

G 2 D(r,<r) = - / dyy 

71 JO 



y 2 + 4 



Mi 



Mi 



(A.l) 
(A.2) 



where Jo(x) is the Bessel function of the first kind. 

dyy 

TT A JO 



G m (r,a) = -^ dyy*—^ 



exp 



-Ml 



(A.3) 



where jo(x) = sin(x)/x is the spherical Bessel function of the first kind. For the special 
case r = we have with z = y/2a/t; 



Gi D (0, a) = - exp(^ 2 ) erfc(z) 
G 2B (0,a) = -^exp(z 2 )Ei(-z 2 ) 



G 3D (0,<r) 



TT 



y/nz 



— exp(z 2 ) erfc(z) 



(A.4) 
(A.5) 
(A.6) 



where erfc(x) is the complementary error function and Ei(x) is the exponential 
integral [56]. 



Appendix B. Derivation of the GME 

The Hamiltonian for the single impurity and the BEC is 

ho = hu B J2 j hj + hUqtiqbq 



3 q 
h I = -JY / (X l a l )\X 

(id) 



J a j) -I 



where ho is the unperturbed Hamiltonian and hi is treated as a perturbation in the 
strong-coupling regime £ <C 1. 



Transport of strong-coupling polarons in optical lattices 



20 



Starting point of the derivation of the GME is the Liouville- von Neumann 
equation [38] 

^ = £««), (B.i) 

where p(t) is the density matrix of the system expressed in the eigenbasis of h . 
The Liouville operator £ is defined by £ = £ + £>ii with C = —i/h[h , ■} and 
Ci = —\/%[hi, •]. The density matrix p(f) is decomposed as p(t) = Vp(t) + Qp(t), 
where V is the projection operator on the relevant part of p(t) and Q = 1 — V the 
complementary projection operator to V . The time evolution of the relevant part Vpit) 
is governed by the Nakajima-Zwanzig equation [41] 

^-Vp(t) = V CP pit) + fdsVCe^QCPpit - a) + VCe QCt Qp(0) . (B.2) 
at Jo 

Specifically, the projection operator V employed in the derivation of the GME is defined 

by Pp{t) = Pb ®VTi B pit) [38,39]. Here, ps is the density matrix of the condensate in 

thermal equilibrium, Trg is the trace over the condensate degrees of freedom and T> is 

the projection operator on the diagonal part. Thus QA(t) = VTrspit) is the diagonal 

part of the reduced density matrix. We note that the trace over the Bogoliubov phonon 

states in the definition of V introduces irreversibility into the system. 

The Nakajima-Zwanzig equation can be simplified given that VCq = CqV = and 

VLfP = for the definitions of Co, £i and V above. In addition, we assume that 

the initial density matrix of the total system has the form p(0) = ps® £a(0) so that 

Qp{0) = 0, and hence the inhomogeneous term in Eq. (1B.2p vanishes. Taking these 

simplification into account and using Tr B Vp = Qa we find from Eq. (IB.2[) that QA{t) 

evolves according to [38] 

%-Q A (t) = f t dsK(s)g A (t-s) (B.3) 
at Jo 

with the memory kernel 

K(s) = VTr B [C^+^d p B (0)] . (B.4) 

The memory kernel to second order in (, i.e. dropping Cj in the exponent in Eq. flB.4p . 
can be expressed in tetradic form as [39,41] 

«hM= E e-^>/ fc - r |(i,{iV q }'|^U,{iV q }}| 2 

n {iV q },{iV q }' 

X COS [(fijj + huj {Nq} - huj {Nci y) t , (B.5) 
with the partition function of the condensate 

Z B = ]T e -^{^> /fesT , (B.6) 

and where Hflij is the energy difference between the impurity configurations % and j, 
and hu)i N \ is the energy of the phonon configuration {iV q }. 
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The explicit expression for the memory function Wij(s) = K^^{s) for the 
Hamiltonian h + hj has been evaluated in [39,40] based on Eq. (1B.4j) . however, we 
here give an alternative derivation of Wij(s) starting from Eq. ( IB. 51) . The evaluation of 
Kii,jj( s ) m Eq. (1B.5P can be separated into a phonon part and an impurity part, where 
the latter is given by 

J 2 ^|(i|4a fc±1 ij)| 2 e- in -* = J% i±l e^. (B.7) 

k 

Thus, for the phonon part we only have to consider operators of the form XjXj±x, which 
can be written in terms of displacement operators D((3) = exp(/%t — (3*b) as 

X}Jt m = l[D(P j!Cl )^, (B.8) 
q 

with /3j iq = M* — M* ±1 q and q the corresponding phase. This allows us to treat 
each phonon mode in Eq. (1B.5P separately, and the problem reduces to the summation 

oo oo 

j2 e- nhuJ ^ kBT \(m\D(P J<(l )\n)\ 2 e [ ^ t(m - n) , (B.9) 

n=0 m=0 

where m and n are phonon occupation numbers of the mode o; q . The matrix elements 
of the displacement operator in the Fock basis are [52] 

(m\D(P)\n) = ^M(3 m - n e-^ 2 / 2 L^- n {\(3\ 2 ) for m > n 

™L (B.10) 
(m\D(/3)\n) = y^(-/3*)"- m e^l 2 / 2 L"- m (|/3| 2 ) for n > m, 

where L^(a;) are generalized Laguerre polynomials. At this point we introduce the new 
variables x = |/3j, q | 2 , y = e 1Wqt , z = Q~ hul i/ k B T ? anc i / = m _ n _ p or fa e case m > n, or 
equivalently I > 0, we find after the substitution of m = n + 1 that the sum in Eq. (1B.9I) 
becomes 

Z=0 n=0 l n + 

To evaluate the sum over n we use the fact that the following relation for generalized 
Laguerre polynomials holds [57] 

to r(n + 7 + l) " 1-z 6XP l Z l- zj^y l \-z) ' 

provided that |z| < 1. Here, is the Gamma function and I 7 (x) are modified Bessel 
functions of the first kind. Using the relation in Eq. (IB. 121) we find that the sum over n 
yields 

T^M-B)^(f^). (B,3, 

For the case n > m, or equivalently I < 0, we have to express L r ^ m (x) in terms of 
~L™~ n (x) in order to exploit Eq. (IB. 121) . Using the relation 

L-'(x) = -^-L' r _ 8 (x) , (B.14) 
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with (x) n = x(x + l)(x + 2) • • • {x + n — 1) the Pochhammer symbol, we find that 

2 j,2(m—n) 2 
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and since 



L-(— \x) 



~n^} (m—n) 



-m, 



(m—n) 



l: 



[(-m)(-m + 1) •••(-«- l)f 



I = m — n m expression (1B.9j) yields 



i=0 n=0 I" ^ V- 



(B.15) 
m!/n!) 2 substituting 

(B.16) 



which is identical to expression flB.111) except for the upper limit in the sum over I. 
Using relation flB.121) again and discarding the double counting of / = we find that 
the total sum in Eq. (IB. 91) is given by 



1-z 



exp 



-x 1 



2z 
1 - z 



To evaluate the sum over / we use the identity [56] 

oo 

h(x)t l = exp 



l=— oo 



X ( 1 
- [t+ - 

2 V t 



and find that Eq. flB.17j) equals 



1 



exp 

1 — z 

which can be written as 



-x 1 + 



2z 
1 - z 



exp 



Xy/z l y v ; 



1 - z \ v z y 



(B.18) 

(B.19) 
(B.20) 

1 _ e -nojq/k B T\-l _ 



Z q exp [ - |^, q | 2 {(iV q + 1)(1 - e^') + iV q (l - e"^*)} 

with iV q = z/(l - z) = (e fo o/ fcflr - l)" 1 and Z q = 1/(1 - z) = ( 
Taking the impurity part and the product of all phonon modes q into account we find 
the complete memory function 
2 r 

exp{- ]T \Pj 



J 



W iJ (s)=25 jii±1 [-\ Re 



q^O 



exp(±io;st) 



(B.21) 



x [(iy q + l)(l-e^ s ) + iV q (l-e- 
where we used Z B = n q Z q . 

Appendix C. Self-trapping 



Impurities immersed in a BEC get self-trapped for sufficiently strong impurity-BEC 
interactions, even in the absence of an additional trapping potential [22-24]. Based on 
the results for static impurities in Section [3] we now show that for the parameter regime 
considered in this paper self-trapping effects can be neglected. 
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As pointed out by Gross [58] the coupled equations describing the impurity and the 
condensate in the Hartree approximation are given by 

^V(r) = -^V 2 ^(r) + (7^(r)|V(r) + K | X (r)| 2 ^(r) (C.l) 

s X (r) = - -^V 2 X (r) + *#(r)| 2 X (r) , (C.2) 

where x( r ) is the wavefunction of the impurity, m a is the impurity mass and e the 
impurity energy To determine whether the impurity localizes for given experimental 
parameters one has, in principle, to solve the coupled equations for ip{v) and x( r ) [24]. 
Alternatively, as suggested in [23], we use the Gaussian mode-function Xj>( r ) in Eq. (jlj) 
as a variational wavefunction for the impurity, with the spread o as a free parameter, 
and minimize the total energy of the system in the regime a C 1, where the linearization 
of the GP equation is valid. For a homogenous condensate, the potential energy of the 
impurity and the BEC is (—£ p ), and thus adding the kinetic energy of the impurity 
yields the total energy 

S{a) = e^a) - £ P {a) = - ^(0, a) . (C.3) 

The impurity localizes if £(cr) has a minimum for a finite value of a, which depends on 
the dimensionless quantity 

/ ^ Wla £p /m a\ 

a = a ~ — — . (G.4) 

g m b £ kin 

In one dimension, we find by asymptotically expanding £(cr) in the limit cr/£ ^> 1 that 
there exists a self-trapping solution for arbitrarily small a' and that 

a lD = V2^- r (C.5) 
a' 

For the two-dimensional case, asymptotically expanding £(cr) in the limit cr/£ ^> 1 yields 
a critical value a' c = 2n, above which self-trapping occurs. The corresponding spread 
of the self-trapping solution o"2D diverges close to a' c , which validates the asymptotic 
expansion of £ (a). In three dimensions, numerical minimization of £{cr) shows that the 
critical value is a! c ~ 31.7, and the corresponding self-trapped state is highly localized 
with a 3 D ~ 0.87^. 

According to Eq. flC.5j) the spread of the self-trapping solution exceeds several 
lattice spacings for the parameter regime considered in this paper. In other words, the 
spread o\d is much larger than the harmonic oscillator length ae in practice, and hence 
self-trapping effects are indeed small. 
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